R Notebook

library(stringr)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.4     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ forcats 0.5.1
✓ readr   1.4.0     
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x readr::col_factor() masks scales::col_factor()
x purrr::discard()    masks scales::discard()
x tidyr::expand()     masks Matrix::expand()
x dplyr::filter()     masks stats::filter()
x dplyr::lag()        masks stats::lag()
x tidyr::pack()       masks Matrix::pack()
x tidyr::unpack()     masks Matrix::unpack()

response <- try(system('~/google-cloud-sdk/bin/gcloud projects list --quiet', intern = T))
projectid <- strsplit(response[2], " ")[[1]][1]

options(na.action = "na.fail") 

source("./dredge_functions.R")

dredge_and_subset <- function(data) {
  model <- lm(
    presence_ratio ~ foraging_niche + trophic_niche + is_noctural + pc1 + pc2 + pc3 + pc4,
    data=data
  )
  dredge_result <- dredge(model)
  summary(model.avg(dredge_result, subset = delta < 10))
}

load_species_data <- function(poolname) {
  filename <- str_replace('species_metrics_##POOL_NAME##.csv', '##POOL_NAME##', poolname)
  
  if (!file.exists(filename)) {
    column_name = str_replace_all('present_in_##POOL_NAME##_pool', '##POOL_NAME##', poolname)
    if (poolname == 'both') {
      column_name = 'present_in_both_pools'
    }
    
    sql <- str_replace_all("
         WITH species AS (
              SELECT 
                  scientific_name,
                  count(*) AS regional_pool_count,
                  countif(present_in_city = true) AS city_count
              FROM model.regional_species_filtered
              JOIN model.taxonomy USING (scientific_name)
              WHERE ##COL_SELECT## = true
              GROUP BY scientific_name
          )
          SELECT 
              species.*,
              body_morphspace.pc1,
              body_morphspace.pc2,
              body_morphspace.pc3,
              body_morphspace.pc4,
              trophic_niche,
              foraging_niche,
              is_noctural,
              ttc.cluster,
              ttc.cluster_half,
              ttc.cluster_double
          FROM species
          JOIN model.taxonomy USING (scientific_name)
          JOIN datasci.taxonomic_trait_clusters ttc USING (scientific_name)", '##COL_SELECT##', column_name)
  
    tb <- bq_project_query(projectid, sql)

    data <- bq_table_download(tb)
    write_csv(data, filename)
  }
  
  data <- read_csv(filename)
  
  data[is.na(data$foraging_niche),]$foraging_niche <- 'Omnivore'
  
  data$foraging_niche = as.factor(data$foraging_niche)
  data$trophic_niche = as.factor(data$trophic_niche)
  data$is_noctural = as.factor(data$is_noctural)
  data$presence_ratio = data$city_count / data$regional_pool_count
  
  data
}

merlin_data <- load_species_data('merlin')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

merlin_data

merlin_result <- dredge_and_subset(merlin_data)

Fixed term is "(Intercept)"

merlin_summ <- model_summary('merlin', 'species', merlin_result)
merlin_summ

birdlife_data <- load_species_data('birdlife')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

birdlife_result <- dredge_and_subset(birdlife_data)

Fixed term is "(Intercept)"

birdlife_summ <- model_summary('birdlife', 'species', birdlife_result)
birdlife_summ

both_data <- load_species_data('both')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

both_result <- dredge_and_subset(both_data)

Fixed term is "(Intercept)"

both_summ <- model_summary('both', 'species', both_result)
both_summ

either_data <- load_species_data('either')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

either_result <- dredge_and_subset(either_data)

Fixed term is "(Intercept)"

either_summ <- model_summary('either', 'species', either_result)
either_summ

Full result

all_species_results <- full_join(full_join(merlin_summ, birdlife_summ), full_join(both_summ, either_summ))

Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")

write_csv(all_species_results, "species_result.csv")

library(ggpubr)

jpeg(width = 1024, height = 1024)
myplot
dev.off()

null device 
          1